Earthquake is a common phenomenon undoubtedly known to humans from earliest times. With the establishment of the plate tectonics theory, scientists started to understand the cause and nature of earthquakes. However, statistical studies appeared much earlier, urged by the necessity of predicting earthquakes.

Gutenberg and Richter realized that the energy released during the earthquake is increased exponentially with the earthquake magnitude.^[1–3] This is the Gutenberg–Richter law. Larger earthquakes occur less frequently, the relationship being exponential (https://earthquake.usgs.gov/earthquakes). In the United Kingdom, for example, it has been calculated that the average recurrences are: an earthquake of 3.7–4.6 every year, an earthquake of 4.7–5.5 every 10 years, and an earthquake of 5.6 or larger every 100 years from (www.quakes.bgs.ac.uk/hazard).

Earthquakes are probably the most relevant paradigm of self-organized criticality (SOC) that can be observed by humans on the earth.^[4–6] The SOC concept was introduced as a possible explanation for the widespread occurrence in nature of long-range correlations in space and time by Bak, Tang, and Wiesenfeld.^[4,7]

In their sandpile model both the random, slow addition of “blocks” on a two-dimensional (2D) lattice and a simple, local, and conservative rule drive the system into a critical state where power law distributed avalanches maintain a steady regime far from equilibrium.

Olami, Feder, and Christensen (OFC) made an important contribution to the SOC ideas by mapping the Burridge–Knopoff spring-block model^[8] into a nonconservative cellular automaton,^[9,10] simulating the earthquake’s behavior. The model results in a power law distribution of avalanches similar to the Gutenberg–Richter law and also reproduces other characteristics of real earthquakes.^[11]

In order to improve resemblance with the geological structure of the earth, and the presence of power law distribution of avalanches, the topology of connections between dynamical units plays an important role. In the literature, OFC models on different topologies have been explored previously. The anisotropic version,^[12–14] random-neighbor version,^[15–18] quenched random graph,^[19,20] small-world,^[21,22] scale-free,^[23–26] random spatial network,^[50] and other different topology models^[50,50,50] have been investigated.

Earthquakes as SOC phenomena, in addition to the lack of success in predictability, have developed the idea that the crust is at a critical state where a minor perturbation can trigger an earthquake of any size and duration, making them inherently unpredictable.^[50] However, foreshocks, aftershocks, and clustering properties indicate the existence of correlation between different events.^{[50,50,50,50,50]} Many seismologists believe that large earthquakes are quasi-periodic,^{[50,50,50,50]} but periodic behavior has appeared in theoretical models only as a special or as a trivial solution,^[50,50,50] or as a result of a phase locking due to periodic boundary conditions (BCs),^[50,50,50] or synchronized regions^[50] in cellular automata.

However, a nontrivial quasi-periodic behavior in the avalanche time series with a period proportional to the degree of dissipation of the system has been discovered by Ramos et al. in 2006.^[50] The correlation of time series of earthquakes can indicate impending catastrophic events and is considered to be possible to achieve prediction.^[50,50,50]

In the aim of improving resemblance with the mechanical model they have introduced two variations in the OFC model. (i) Thresholds are distributed randomly following a narrow Gaussian distribution of standard deviation. (ii) Instead of assuming infinitely accurate tuning, a constant and finite force continually drives the system, but keeping the separation of time scales (relaxations are considered to be instantaneous) which is more realistic.

It is natural to think about whether there are quasi-periodic phenomena in other OFC models. And trying to find the mechanism for the emergence of quasi-periodicity.

In the present work, we first investigate the quasi-periodic behavior in a random spatial network. It is not like the original OFC model which exhibits SOC behavior for a wide range of dissipation values and the power-law exponent depends on the dissipation value, SOC can be excited only in the approximate conservative or the network is denser case.^[9,50] There is not SOC on sparse spatial networks, due to the fact that spatial networks include many modules and the modular structure hinders the spreading of avalanches in the whole network. Connection degree is increased which enables energy transfer over a long range. The long-range energy transfer overcomes the effect of local modularity and SOC can be reached. An intuitive scenario is that there are no quasi-periodic phenomena in nonconservative random spatial network.

However, by analyzing the autocorrelation function and the fast Fourier transform for the avalanche time series, we show that there is a quasi-periodic phenomenon in the random spatial network. And we obtain a simple rule for a period that is proportional to the choice of unit time and the dissipation of the system. Although the periodicity of the avalanche does not depend on the criticality of the system or the average degree of the system or the size of the system, there is evidence that it depends on the time series of the average force of the system.

We have also explored other models in the present study. Cycle graphs and two-dimensional lattice models are in good agreement with the above laws, but small-world models, scale-free models, and random rule graphs do not have periodic phenomena.

The random spatial earthquake model was introduced by replacing the lattice model, in order to improve randomness and locality, which is more realistic.^[50] In which sites are randomly placed on a plane and are connected locally. The network is built by two simple rules: (i) Randomly pick N sites on a square, whose width is L; (ii) Two sites are connected if the distance between each other is less than the connection radius r_c.

The system average degree is , where k_i is the degree of the i-th site, that is, the i-th site is linked to k_i neighbors. For the given values of L and N, the average degree can be changed by tuning the connection radius r_c.

Initialize all sites to a random value between 0 and F_th. When a site F_i reaches the threshold F_th, and a fraction of its force is redistributed to its neighbors, 1 where “nn” stands for the collection of nearest neighbors to node i. The parameter β ( ) controls the level of dissipation of the system and takes values between 0 and 1. If β=1.0 or α=1/4 the system is conservative. According to the rule (1), the force of an unstable site at the boundary is averagely distributed to each of its nearest neighbors. Therefore, we use free boundary conditions in this model.

The toppling of one site triggers an avalanche, that is, neighbors of this site may become unstable and toppling propagates in the network. The avalanche is over until all of the sites are below F_th. The number of toppling sites during an earthquake is defined as the earthquake size S. Then the driving to all sites recovers, a new avalanche is triggered. In a system of SOC, the distribution of earthquake sizes is a power-law function. We still adopt the above two changes: i) Thresholds are distributed randomly following a narrow Gaussian distribution of standard deviation σ=0.001. ii) A constant and finite force continually drives the system, . The modified model is more similar to the mechanical model and more convenient to calculate time series.

We use the natural unit time for simplicity, that is an isolated block would need 10⁴ steps from zero to reach the threshold for , the time is .

The autocorrelation function is defined as 2 where corresponds to the avalanche time series, is the avalanche size average.

The original OFC model exhibits SOC behavior for a wide range of dissipation values and the power-law exponent depends on the dissipation value, and the original random spatial earthquake model exhibits SOC only in the approximate conservative or the network is denser case.^[9,50] Figure 11 shows the avalanche size distribution in finite velocity for the OFC model (see Fig. 1(a) and Fig. 1(b)) and the random spatial model (see Fig. 1(c) and Fig. 1(d)). The model that makes this change is called the finite velocity model, in which a constant and finite force continually drives the system, and keeping the separation of time scales (relaxations are considered to be instantaneous).

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. Distribution of earthquake size for the finite velocity OFC model (see Fig. 1(a) and Fig. 1(b)) and the finite velocity random spatial model (see Fig. 1(c) and Fig. 1(d)). Figure 1(a) and 1(b) show the avalanche size distribution for different system size N and different dissipation α on the finite velocity OFC earthquake model with α=0.2 and N=1282, respectively. Figure 1(c) and 1(d) show the avalanche size distribution for different system size N and different dissipation β on the finite velocity random spatial earthquake model with β=0.8 and N=1282, respectively. The power law behavior of the avalanche distributions is not more robust in this finite velocity random spatial model than in the finite velocity OFC one.

Figure 1(a) and 1(b) show the avalanche size distribution for different system size N and different dissipation α on the finite velocity OFC earthquake model with α=0.2 and N=128², respectively. Figure 1(c) and 1(d) show the avalanche size distribution for different system size N and different dissipation β on the finite velocity random spatial earthquake model with β=0.8 and N=128², respectively. The largest avalanche size in the OFC model is about 10⁴, which is close to the system size N. In the random spatial models, the distribution of avalanche size depends on the dissipation parameter β, but the distributions obviously diverge from a power-law function. Although the system size is N = 128², the largest size of avalanche is very small. The power law behavior of the avalanche distributions is not more robust in this finite velocity random spatial model than in the finite velocity OFC one.

Next, we explore whether quasi-periodic phenomena will occur in random spatial earthquake models by calculating correlation functions. The system has parameters and N=128². The analysis of the autocorrelation function (see Fig. 2(a)) displays a strong correlation between avalanches. Some of the peaks for σ=0.001 are shown in Fig. 2(a). The maximum height of the peak is less than 0.01. Although the height of these peaks is relatively small, the position of these peaks indicates that the system has a quasi-periodic behavior with a period proportional to the degree of dissipation and the unit time 3

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. The quasi-periodic behavior on different earthquake models. The first row show the analysis of the autocorrelation function for the avalanche time series of (a) the random spatial earthquake model with , (b) the OFC earthquake model, and (c) the earthquake model on cycle with k = 4. Each type of graph or network is simulated for the system size N=1282. Different curves correspond to different degrees of dissipation or α. The position of the peaks indicates that the system has a quasi-periodic behavior with a period proportional to the degree of dissipation: . For visual clarity, curves α=0.10 and α=0.15 have been shifted along the x axis in the earthquake model on cycle, and , respectively. The second row shows the fast Fourier transform (FFT) for the avalanche time series of (a) the random spatial earthquake model with , β=0.6 and N=1282 ( ), (b) the OFC model with α=0.15 and N=1282, and (c) the earthquake model on cycle with α=0.15, N=1282, and k = 4. The quasi-periodic behavior is robust T=1/f and the period is consistent with that of the autocorrelation function.

Surprisingly, the period under this random spatial earthquake model is consistent with that of the OFC model (see Fig. 2(a) and Fig. 2(b)). Although the height of the peak is also small, the position of the peak shows that the quasi-periodic behavior is obvious. However, the curve of OFC model is higher and smoother than that of random spatial earthquake model, because the random spatial system has more perturbations. The randomness and locality of the random spatial network make the noise of the correlation function more stronger than that of the 2D spring-block model.

We also see that the peak value of the correlation function decreases first and then increases with the increase of the dissipation coefficient in the random spatial earthquake model (see Fig. 2(a)). However, the peak value of correlation function increases and then decreases in the lattice model (see Fig. 2(b)). The height of the peaks shows a non-monotonous variation with the dissipation degree not being well understood yet.

It is not like that periodicity is as robust as criticality in OFC model, periodicity is independent of criticality in this model. The original OFC model exhibits SOC behavior for a wide range of dissipation values and the power-law exponent depends on the dissipation value, SOC can be excited only when the average degree is large and the system tends to conservation in this random spatial model. From the above research, we can see that this random spatial model, like the original model, has quasi-periodic behaviors under various conservative parameters. Therefore, periodicity and criticality are not dependent on each other.

Besides, the quasi-periodic behavior is robust for different system average degree and different system size N on the random spatial earthquake model. Nevertheless, the peak is getting lower and lower as the system average degree increases (see Fig. 3(a)). Although the position of the peak is just slightly offset (see Fig. 3(b)), the quasi-periodic behavior is also robust for different system size N.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. The analysis of the autocorrelation function for (a) different system average degree with N=1282 and β=0.60, (b) different system size N with and β=0.60. The fast Fourier transform (FFT) for (c) different system average degree with N=1282 and β=0.60, (d) different system size N with and β=0.60. Except for large errors when N is large, the quasi-periodic behavior is also robust for different system size N or system average degree .

Moreover, we computed the fast Fourier transform (FFT) for the avalanche time series of both the random spatial earthquake model with and the OFC earthquake model as the second row in Fig. 2. Each type of graph or network is simulated for the dissipation degree β=0.60 or α=0.15. The quasi-periodic behavior is robust, (T=1/f, f is the frequency at the peak) and the period is consistent with that of the autocorrelation function. We also computed the fast Fourier transform (FFT) for different system average degree k and different system size N on the random spatial earthquake model in Fig. 3(c) and Fig. 3(d). Except for large errors when N is large, the quasi-periodic behavior is still robust.

We also discussed the mechanism of periodic phenomena. In general, the inherent complexity of SOC on networks makes analytical investigations almost always impossible. Although the periodicity of avalanche does not depend on the criticality of the system or the average degree of the system or the size of the system, there is evidence that it depends on the time series of the average force of the system.

The analysis of the average force of the system , corresponds to the time series (gray lines in Fig. 4), displays a strong correlation between the system average force: 4 Surprisingly, the period is consistent with that of the autocorrelation function and it is depending on the time series of the average force of the system (see table 1). The fitted equation is 5 where a, b, y₀, and x₀ are fit parameters, see table 1 for fitting details (red lines in Fig. 4). Except for large errors when β is small in the random spatial model, the period T is consistent with that of the autocorrelation function.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. The analysis of the average force for the time series of (a) the random spatial earthquake model with and N=1282, (b) the OFC model with N=1282, and (c) the cycle model with k = 4 and N=1282. Different curves correspond to β=0.40 (α=0.10), 0.60 (α=0.15), and 0.80 (α=0.20) from up to bottom. Gray lines represent the average force for the time series and red lines are the fitting of the average force data of the system (see table 1 for fitting details).

Table 1.

Table 1.

Table 1. Fitting details corresponding to Fig. 4. The fitting equation is Eq. (5). . Models β (or α) y0 x0 b a T Random spatial model 0.4 0.50022 17645.28997 1475.5896 0.00185 2951.179 Random spatial model 0.6 0.50065 8060.66579 992.4178 0.000487 1984.836 Random spatial model 0.8 0.50147 37670.52985 964.8767 0.0009638 1929.753 Lattice model 0.1 0.50164 7243.2208 2988.81158 0.00602 5977.623 Lattice model 0.15 0.50954 614.73525 2003.725 0.01609 4007.45 Lattice model 0.2 0.5512 –48899.38938 1142.83261 0.02087 2285.665 Cycle model 0.1 0.50004 7902.04672 2987.4371 0.01366 5974.874 Cycle model 0.15 0.50008 1691.95799 2001.83301 0.0171 4003.666 Cycle model 0.2 0.50007 –242.57117 999.64596 0.0023 1999.292

Table 1. Fitting details corresponding to Fig. 4. The fitting equation is Eq. (5). .

Moreover, we find that the rule also holds for cycle graphs with k = 4 by analyzing the autocorrelation function for the avalanche time series in Fig. 2(c) and the average force for the time series in Fig. 4(c). On cycles, individuals interact with their k/2 nearest neighbors on either side. Thresholds are distributed randomly following a Gaussian distribution of standard deviation σ=0. Instead of showing non-monotonic variations in peak height, peaks almost always reach a maximum of 1.0. The peaks for small σ are extremely narrow. It is almost coincident covering for different dissipation parameters. For visual clarity, curves α=0.10 and α=0.15 have been shifted along the x axis in the earthquake model on cycle, and , respectively.

To illustrate the period depends on the time series of the average force of the system, we also analyze the average force time series of the anisotropic model, random rule graph, small-world network, and scale-free network, and find that these models do not exhibit quasi-periodic behavior. This is consistent with that we calculate the autocorrelation function to the avalanche time series for these models. Presumably because they have a larger variance of the degree distribution.

These numerical results indicate that periodicity is determined by the system average force. If there is quasi-periodic behavior in the time series of system average force, then there is also quasi-periodic behavior in the time series of avalanche. And the period is proportional to the choice of unit time and the degree of dissipation.

We have investigated how structure affects the quasi-periodic behavior of earthquake model. More specifically, we have shown that the quasi-periodic behavior is robust for structured models by the analysis of the autocorrelation function and the fast Fourier transform, including random spatial models, 2D lattice models, and cycles. We also obtained a simple rule for a period that is proportional to the choice of unit time and the dissipation of the system. It is not like that periodicity is as robust as criticality in OFC model, while the periodicity is independent of the criticality in this model. Next, we have explored other models. Although random spatial models, cycle graphs, and 2D spring-block models are in good agreement with the above laws, the small-world models, scale-free models, and random rule graphs do not have periodic phenomena. We have also explored the mechanism by which periodic phenomena occur. The periodicity of avalanche does not depend on the criticality of the system or the average degree of the system or the size of the system, however, there is evidence that it depends on the time series of the average force of the system. Since the quasi-periodic behavior exists in the time series of the system average force, the quasi-periodic behavior also exists in the time series of avalanche, and the value of the period is consistent. If the periodicity and self-organized criticality of earthquakes are taken as the basis, the random spatial model and the two-dimensional lattice model are most suitable.

We would like to thank Lin Zhou for helpful discussions.